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dependence of the stability range upon Igiz, where v is the viscosity. This effect depends upon 
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is brought to zero. These Rndings show how viscous effects play a role in determining inertial 
properties of shell models and give some hints for understanding the effects of viscous dissipation 
upon real turbulence. 
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1 Introduction 


Kolmogorov’s classic 1941 paper [1] asserts that in the limit of large Reynolds number, the inertial 
range behavior is independent of the form and magnitude of the dissipation. Many approaches 
such as large eddy simulations (LES, [2]) or simulations with hyperviscosity [3] have been based 
upon this idea that the viscous subrange (VSR) simply serves as a drain to remove the energy 
carried down through the inertial subrange (ISR). However, recently Leveque and She [4] have cast 
doubt upon the independence of the ISR on the VSR energy sink. They suggested that energy 
cascaded through the inertial range tended to be reflected at the dissipative scale. This reflected 
current could reappear in the integral range and thus could affect the details of integral range 
behavior. In particular, through this mechanism ISR scaling exponents could become Reynolds 
number dependent as also suggested by the log similarity model [5]. 

A variety of approximate models have been developed in recent years to test our ideas about 
turbulence. These models are shell models in which the velocities are put on a fractal series of 
shells in wave vector space. The different shells are characteristized by wave vectors of the form 

kn = kor ( 1 ) 

with A being the ratio between shells. One of these shell models is the so-called GOY model, 
introduced by Gledzer [6], Yamada, and Ohkitani [7], which involves a forcing at large scales, 
a cascade of energy through a set of shells in an inertial range, and a dissipative dominance of 
shells beyond the inertial range. (For some more recent studies see for example [8, 9].) The GOY 
model is perfectly set up to test the hypothesis that inertial range behavior is independent of what 
happens in the VSR. In contrast to full numerical simulations [10] and to reduced wave vector set 
approximations (REWA) [11], even an semi-analytic approach will be possible. 

In this paper, we test Kolmogorov’s basic hypothesis - the independence of ISR quantities on 
the VSR - by looking at the static properties of GOY model behavior. As Rrst pointed out by 
Biferale et. al. [12] this model has a static solution for certain values of parameters. In this paper, 
we find the natural static solution (section 2), define the dependence of this solution upon the 
model parameters (section 3,4) and do linear stability analysis to determine the range of stability 
of the static solution (section 5). 

The major conclusion of this paper is that Leveque and She [4] were right for the GOY model 
at least. The inertial range behavior of the static solution does depend upon the strength and 
structure of the dissipation. As we develop below, both the static solution and its range of stability 
are - for small values of the viscosity - periodic functions of the logarithm of viscosity, with the 
period being proportional to the logarithm of A. The fundamental source of this periodic behavior 
is the meshing between the onset of dissipation and the shell structure. Thus, for example, the 
dissipation of energy has a different value when the dissipative term becomes first large at a given 
shell or if conversely it becomes large between shells. As the range of stability depends on the 
viscosity and as inertial range scaling corrections to K41 (commonly called 6(p = (p — p/3 for 
the p-th moment of the velocity) set in smoothly at the borderline of stability [12], this suggest 
that also the scaling corrections 6(p depend on viscosity, maybe even in the fully chaotic state, as 
already discussed in [4, 9]. Since the shells are unphysical, this mechanism for viscosity dependence 
is also unphysical. In fact, this periodic behavior disappears when the shell ratio goes to unity. 
However, some roughly analogous mechanism for coupling between inertial range and dissipative 
range behavior might even exist in real turbulence and thus contradict some portion of the lore 
surrounding the K41 paper. This item will be discussed in our conclusions (section 6). 


2 The GOY model and its static solution 

The basic idea of the GOY model is to represent the wavevector space by N geometrically scaling 
wavevectors n = 1, 2, . . ., V as in equation (1). A complex number Un represents the typical 
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velocity on scale k~^. It is coupled to next and next nearest neighbors. Large scales are forced. 
Viscous damping becomes effective on small length scales. To be more precise, the GOY toy 
dynamics for the three-dimensional Navier-Stokes equations is given by 

jUn = i [Cn{U)]* - l^kP„Un + (2) 

Here the cascade term, Cn(U), is a bilinear expression in U which links neighboring shells. In the 
GOY model it has the form: 

On(G) = knUn-\-lUn-\-2 — 1 — 1 Gn -|-1 (1 €)kn — 2Un — lUn — 2 (3) 

while the terms in F and v are respectively representations of the large-scale forcing and the 
small-scale viscous dissipation. The parameters in the model are the scaling factor A between the 
shells, the parameter e determining the ratio among the three cascade terms, and the viscosity v, 
which can be considered as inverse Reynolds number. Most of the calculations in this paper are 
carried out for “normal” viscosity {p = 2). Sometimes we examine the effect of hyperviscosity by 
looking at p > 2. Normal viscosity is to be assumed unless stated otherwise. The forcing F and 
the wavenumber ko determine length and time scales, we make the conventional choices ko = 
and T = 5 • (1 -|- i) • 10“^. In the inviscous, unforced case, equation (2) conserves the total energy 
Xln |GnP/2. A second conserved quantity is |i7„p(e—1)“" which can be identihed with helicity 
if A = 1/(1 — e) [9]. The number of shells N is chosen big enough so that the energy on the last 
shells is practically zero. Much work has been done in recent years to understand the dynamics of 
equation (2) [ 8 , 9, 13]. 

For the static case the complex equation (2) can be set in real form. A simple change in the 
phase of the Un can eliminate all the i’s in the equation and leave the result real. There is a 
periodic behavior in n built into the inertial range, in which the behavior on the different shells 
repeat with a period three [9]. The change in phase is done with the basic period three [13] in the 
form Un = UncF’'^ ^ F = /e®^. By taking $ = 7 r /4 we make / real and equal to 5 • \/2 ■ 10“^ in 
the standard case. From dynamical calculations we Rnd that the static solution of the complex 
dynamical equation ( 2 ) picks the phases 

( 97 r /8 for n = 3, 6 , 9,... 

</n = < 7 r /4 for n = 1,4, 7, . .. (4) 

[ 7r/8 for n = 2, 5, 8 ,... 

So we also choose this very phases and thus equation (2) is replaced by the real equation 

^Un = -Cn{u) - vkP^Un F fbn,A- (5) 

with Cn(u) defined in equation (3). Equation (5) is by the way Gledzer’s original model for 
two-dimensional turbulence [ 6 ]. 

A phase choice like (4) can turn every static solution into real form. This is an exact result for 
nonvanishing viscosity [14]. In the examined cases, for the above chosen phases, all are positive, 
apart from M 3 which is negative and very small in modulus, (us goes to zero with 12 0 .) 

In the static case, the left hand side of eq. (5) is set to zero. The result is a fourth order 
difference equation which will then require four pieces of boundary data. Our solution is obtained 
by using the boundary conditions mq = m_i = mjv-i-i = 'Un +2 = 0. Other boundary conditions are 
possible, e.g., one could think of fixing mq = 1 rather than forcing the system, but qualitatively 
they lead to the same conclusions as those boundary conditions chosen here. 

In order to Rnd the solutions m„ of these N real equations, we Rrst get a rough solution by 
iterating forward the time dependent equations (2) for some small value of e. After a while, the 
velocities approach a reasonably constant set of values. This roughly static solution is then reRned 
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by Newton’s method. After the exact static solution is obtained for some set of parameters, (i', X, e), 
then we continuously vary one or the other parameter(s) in small steps to obtain neighboring 
solutions. Once solutions are obtained over a reasonably large range, we can use the theory 
outlined below to extend the solution to smaller values of the viscosity and to arbitrarily large 
numbers of shells. Once there are enough shells, adding more just extends the VSR but does not 
change the values of the velocities in the other shells. 

Since these static equations are non-linear and of high order, one expects many solutions. In 
this paper, we focus upon the particular solutions which evolve naturally after a transient time. 
They all have the same character as those in the base case of A equal to 2 and small values of e. 
These solutions exist over a wide range of parameters but do not exist everywhere. For example, 
Rgure 1 shows the region in A and e for which one can find solutions to eq. (5) for the standard 
case ly = 10“^, / = 5 • \/2 ■ 10“^. The number of shells N is picked sufhciently large. This figure 
also shows the stability of our solution, as indicated by linear stability analysis. (See our work in 
section 5 below.) The lines on the figure define the places where the first stability eigenvalues pass 
from their stable to their unstable regions. Thus, the regions to left of these lines are the stability 
domains for the static solutions. The left solid line holds for the stability of eq. (2), the middle for 
that of eq. (5). In general, instabilities in the phase (j)n occur for smaller e than for the modulus 
Un, SO the former is most frequently left to the latter. 

Figure 1 is quite complex, describing many processes at once. We summarize this figure by 
describing what happens along the cut on the figure at A = 2. For small e there is a stable static 
solution of eq. (2). This solution continues to exist all the way up to Cgnb ~ 0.44990622582 where 
there is a saddle node bifurcation [15]. At this bifurcation point one eigenvalue (in this case the 
largest) of the Jacobian of equation (5) becomes zero and a unique curve of fixed points passes 
through the bifurcation point which lies on the side e < Cgnb [15]. The static solution of (5) is 
stable over the entire range ]0,es„j]. Yet this is not the general case. For different A figure 1 
shows that there are many cases where equation (5) becomes unstable as well. For A = 2 equation 
(2) undergoes a Poincare-Andronov-Hopf (short: Hopf) bifurcation [15] at e Ri 0.36987669663 Ri 
0.3699. This transition was analyzed by Biferale et al. [12]. The next most unstable mode turns 
positive at Ri 0.39560317884 Ri 0.3956. If one requires as did Gledzer [ 6 ] that the solution be real, 
then these instabilities cannot appear. Thus we see some of the major process which will affect 
the transition to chaos in the GOY model. 

There is another striking feature of figure 1. It contains many wiggles. Even more pronounced 
wiggles are seen in £stab as a function of Igiz, see figure 12 below. Wiggles of this kind were first 
observed in a numerical study carried out by one of us [14]. One of the primary purposes of this 
paper is to explain these oscillations. As we shall see below they are a consequence of the shell 
structure [16]. 

Next we describe the qualitative nature of the solutions obtained for very high values of N and 
small viscosity. We follow the classical analysis and divide our description into three parts: 

a. A forcing range, n = 1, 2, 3, 4. 

b. An inertial range. In this ISR, the viscous term is much smaller than any one of the different 
cascade terms. These cascade terms then cancel against one another. 

c. A dissipative range. In the VSR, the viscous term cancels again the single largest cascade 
term-and the other two cascade terms are much smaller. 

Start with the inertial range. Here one can obtain [9] a simple asymptotic solution for the 
product: 

Sn — kn — lUn — lUnUn-\-l (6) 

for n > 4 of the form 

S„=A + B(e-ir. (7) 

Here A and B determine the fluxes of the two conserved quantities. We look at the case in which 
the second term in equation (7) decays rapidly with n. In this case, the energy flux 

Jfi — ^n'^n'^n-t-l'^n-t-2 4“ (1 <^)^n — I'^n — 1'^n'^n-t-1 (8) 
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gives the crucial transfer of energy (and information) from one region of n to another, namely 

= (2 - e)A. (9) 


The other flux [9] 

— (^ 1) (^n —I'^n —l'^n'^n + 1 '^n'^n + 1'^n + 2) 

is 

T„ = (2-e)5 ( 11 ) 

and carries no additional information, as, see below, equation (15) Jn/Tn = A/B = (1 — e)^ holds. 

So far we have described two constants A and B of the four constants of integration which 
describe the solution in the ISR. The two additional arbitrary constants arise from the symmetry 
of the ISR equations under a period three multiplicative modihcation of The ISR solution is 
invariant under the multiplicative modihcation of solutions in the form 

i aUn for n = 4, 7, 10, .. . 

l3u„ for n = 5, 8 ,11,... (12) 

jUn for n = 6 , 9,12,... 

Since the modihcations of equation (12) also change A and B, they are not independent of the 
previously dehned constants of integration. However if one demands that the product of the three 
changes be unity: 

1 = a/Sj (13) 

then this freedom will not change the value of T or B. By demanding the constraint of equation 
(13) we then have four independent constants of integration, which we can pick to he A , B, l3 , 
and j. 

Notice the nature of the information transfer through the inertial range. Our boundary con¬ 
ditions hx four parameters. The value of B is of minor importance, as, when calculating as 
a function of the four constants A, B, {3, and 7 , it will not depend on B for large n. However, 
both the value of T - which determines the energy flux and equally the values of /3 or 7 - which 
determine the ratios usk+ 2 /'U 3 k and usk+i/usk - have an entirely non-decaying effect throughout 
the inertial range. Thus these three parameters transfer information without loss up and back over 
the whole ISR. 

To set the four parameters, the system sets two conditions on each boundary of the ISR. To 
understand the integral scale, we note that M 3 is very small and in fact goes to zero with 12 as 
the viscosity approaches zero, as mentioned above. We thus neglect M 3 in all the static equations 
with n greater than 2. The resulting simplihed set of equations is essentially the static version of 
equation (5) 

0 =-Cn{u) - i2kPu„ + f6„^4 (14) 

with modihed integral scale boundary conditions M 3 = M 2 = 0. These boundary conditions are 
similar to those suggested by Biferale et al. [12]. One of these boundary condition sets the ratio 
of A and B. To make equation (7) consistent with M 3 = 0, i.e., 54 = 0 from equation ( 6 ), we must 
take 

B = A/{1 - e)^. (15) 

The other boundary condition for the inertial end is obtained from equation (14) applied to n = 4. 
Since M 3 is negligible, we see 

f = k4U5U6. (16) 

Notice that equation (16) uses the applied forcing to Rx a boundary condition upon the inertial 
range but that it does not in itself fix the value of the energy flux. To know that flux, we would 
need to know M 4 / - which is the energy flux. However, we shall not know M 4 until we have solved 
the entire problem including the boundary condition at the end of the VSR. 
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Deep in the VSR, we demand that the velocity decreases rapidly with shell number. In fact, as 
we shall develop just below, for very large n we can obtain the very rapid decrease of the form: 

lnM„ = —const X (knY (17) 

where g is greater than one. Its actual value is of no importance for our argument. To see how this 
form arises, assume a decrease of this form. Then the (1 — e)-term dominates the cascade terms 
and balances against the viscous terms. We then Rnd, assuming that (1 — e) is positive: 

lnM „_2 + lnM„_i = lnM„ + ln(4i^A;„/(l - e)). (18) 


The last term is asymptotically negligible in comparison to the others. Throw it away, and one 
finds a very simple solution to the resulting Fibonacci equation [17], namely 


lnM„ 


—a 


1 +Vs 

2 


+ h 


1-Vs 


(19) 


where a is positive. Thus, the second order difference equation (18) has a two parameter family of 
solutions which die off rapidly with n. The two parameters a and h are set by fitting onto the inertial 
range solution. Note however, that we have gone from a situation in which we generally have four 
constants of integration to one in which we have but two in the far reaches of the VSR. Where have 
the other two gone? They were, in fact, set by the condition that must decrease rapidly with 
n. This condition eliminated two growing solutions. The demand that the inertial range solution 
must not excite the growing modes in the VSR then sets the two remaining parameters in the ISR 
solution. In this way, all parameters in the ISR solution are set. In particular, the energy flux is 
set by this interplay of ISR and VSR behavior. 


3 Viscosity p and shell distance A dependence of the solution 

We now explore the n and A dependence of the static solution. We saw in the last section that 
the solution was determined by detailed balancing of the various cascade terms and the dissipation 
term on the boundary between the VSR and ISR. The result (and hence the values of the flux) 
will depend in detail upon how the dissipation ’cuts in’. The detailed values of the integration 
constants will be different depending upon whether the dissipation first becomes important on a 
shell for which n equals 3A; or 3A; + 1 or 3A; + 2. The result will return to the same form only after 
the value of n is decreased by a sufficient factor to bring the solution to the same form three shells 
down. Since the cascade term will increase by a factor A over three shells while the dissipation 
term will increase by a factor A® over the same interval, then decreasing by a factor of will 
bring the system back to the same balance of terms three shells deeper down. To see this we plot 
in figure 2 the solution u„ for a succession of cases in which v differs by a factor of A'^. The picture 
shows the same transition from ISR to VSR occurring again and again for the different values of 
n - only for smaller n the transition occurs on higher shells. 

If we vary v by some factor which is not an integer power of A'^, the detailed structure of the 
wiggles in figure 2 will be changed because there is a change in the way the integral scale merges 
into the VSR. To see how this variation works, we look at the standard case A = 2, e = 0.3 and vary 
the viscosity. In figure 3, we show the result of this variation upon the important determinants 
of ISR behavior: the energy current and the ratios usk+ 2 /'U 3 k and usk+i/usk- Each of these 
quantities, for 3k in the ISR, becomes independent of the subscript. Note the periodic oscillation 
of these quantities with Ig As we have already discussed the period is Ig A'^. For the hyperviscous 
case p = 4, the same argument gives a period Ig A^®; in general, it is Ig 

Note the large relative amplitude of the oscillation. For figure 3 we have Ri 0.35. To 

understand this result, we calculate the dependence of on A for fixed n, see figure 4. This 
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oscillation comes from the dissipation of the energy current J which takes place at the high-n edge 
of the ISR within a hxed k-range AA; of, say, about one decade. If there are enough shells in this 
range, i.e. if A is close to 1, then J can precisely adjust to the supplied viscosity. If, on the other 
hand, the shells are sparse, the total dissipation and thus the total flux J sensitively depends on 
the exact positioning of the shells which are situated within the hxed range Ak. We can analyze 
the amplitude of the oscillations when the number of shells in this range is not too small. Since the 
motion caused by a variation in iz can at most move one shell into or out of range, we may expect 
the relative variation in dissipation caused by the motion to be the inverse of the number of shells 
in the range. The number of shells An within Ak is An = lg(AA;)/lg A. Therefore we expect 


—- — = —Ig A 

J An lg(AA;) 


( 20 ) 


When we look at the numerics, we hud order of magnitude agreement with equation (20). For 
A = 2, we may expect an oscillation AJ/J Ri Ig 2 Ri 0.3, as observed. But we also observe some 
non-analyticity as A goes to 1. Such non-analyticities often arise in continuum limit problems, see 
e.g. ref. [18]. 


4 Large Forcing solution 

So far all of our analysis has been performed for large-A systems. As A got close to one we increased 
the number of shells so that, in effect, we were always working in the limit of inhnite shell numbers. 
However, the matching between shells and dissipation which we are discussing in this paper is very 
well illustrated by considering the case of the shell number N hxed. If we then make the forcing 
constant / very “large”, nearly all of the dissipation will occur on the last two shells. This gives 
us the dehnition of the critical / = /g, beyond which we consider / as large. We dehne fc by the 
condition 

nkjfujf ~ fcU4. (21) 

With Kolmogorov scaling mjv ~ ‘^^{ki^/k 4 )~^l^ and M 4 ~ \/TF^ on dimensional grounds we obtain 

/, = fc-]iz2A(8^+4)/3 (22) 

as critical value for /. 

In Rgure 5 we display the / dependence of M 4 , M 5 , and ug for the standard case A = 22, A = 2, 
e = 0.3, iz = 10“^. Since the coefficients A and B have the same /-dependence, we can restrict 
ourselves to theses three all higher of the same triality have the same scaling dependence 
as the corresponding basic velocity M 4 , M 5 , or ug. We observe two quite different scaling regimes 
in figure 5. Indeed, the crossover between them agrees with what we calculate from equation (22), 
namely fc = 2 . 8 . 

For small f < fc, all Un (n 3) scale as Un ~ f^^"^, apart from wiggles of period A® in Ig /. The 
Un ~ f^^"^ scaling follows from dimensional analysis. Lengths are measured in terms of [A;/^], times 
in terms of [fkif\~^l‘^. So velocities are measured in terms of [Ajq ^/^/i/2j ^nd thus scale as ~ f^^"^■ 
From equation (5) it follows that M 3 = vk\U\ju 2 ~ ~ 1. The wiggles on the underlying scaling 

law mirror the above explained fluctuations in v for fixed /. As [v\ = length/^/time = [Ajq 
the period A'^ in Ig iz means a period A® in Ig/ which exactly is what we observe. 

For large f > fc the dissipative mechanism is quite different from the Navier-Stokes case, yet 
it is similar to large eddy simulation-type models (see e.g. [ 2 ]) with some energy drain for large 
wavevectors. We start from K41 scaling, which can be expressed as 

'^n — ^'^n-t-3* (^^) 

Thus, we focus upon predicting the /-scaling of M 3 j,_|_i, M 3 j,_|_ 2 , usk- Surprisingly, it depends on the 
total number N of shells. The reason is that due to the n ^ n + 3 symmetry it matters whether 
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exclusively complete triples (m„, Mn+i, Mn+2) fit between the forced shell n = 4 and the dissipative 
shells n = At, At — 1, or whether the last triple has to be truncated to (m„, Mn+i) or even to (m„). 
Note that the total number of complete triples between n = 4 and the dissipative shells At, At — 1 
does not matter. We thus have to distinguish between the cases At = 3A; + 1, At = 3A; + 2, and 
N = 3k. The constraints on the m„’s are the same in all three cases. Equation (16) provides a 
small-n constraint. The large n-relations are obtained from equation (14) applied for n = N and 
n = At — 1. The resulting constraints are 


Mjv ~ un-iUn-2, (24) 

Mjv-i ~ umum-2, (25) 

but via the n ^ n + 3 symmetry they translate to quite different /-scaling for M 4 , M 5 , and ug for 
the cases N = 3k + l,3k + 2,3k, which are displayed in table 1. 

The spectrum for the case At = 3A; -f 1 = 22 is shown in Rgure 6 . It again reveals the period 3 
structure of the GOY equations. As follows from table 1, the amplitude of this oscillation linearly 
grows with /. 

The analysis of this chapter again reveals how the detailed structure of the VSR and the stirring 
subrange SSR determine the structure of the whole solution throughout. One should not be too 
surprised: Since equation (5), as discussed above, can be considered as a difference equation in 
the variable n. Again its solution is determined by the left and right boundary conditions for the 
equations for n = 1, 2 and n = N, N — 1. 


5 Stability analysis 


Up to now we only discussed the existence of the solutions and its dependences on various param¬ 
eters, but not whether the solution is linearly stable. Yet the case treated in recent publications 
[ 8 , 9] is the unstable one, since it is only in this region that the GOY equations make sense as a 
model for turbulence. We would like to explore how unstability is achieved and how it depends 
on the various parameters. The first stability analysis for the GOY systems was performed by 
Biferale et al. [12] for one set of parameters V = 19, A = 2, = 10“®, ko = 0.05. They find that 

instability in the GOY equations (5) sets in via a Hopf bifurcation (i.e., a complex conjugate pair 

of eigenvalues passes through the imaginary axis) at e = egtah ~ 0.38. We find that this mechanism 
is rather general. 

To discuss the stability of the GOY equations, we return to the consideration of equation (2). 
We are interested in the stability properties of the complex solution Un = Wn exp (i/n)- Variations 
in Un are split in those in the modulus 

Un = + hUn (26) 

and those in the phase 

= + (27) 

The stationary phase (f)^n^ is given by equation (4). The stability matrix equation for the modulus 
is readily obtained from equation (5) as 

— = Yj{Dnm T CnmffUnn. (28) 


Here Dnm and Cnm 
response defined as 


and 


are respectively the dissipation and the cascade contributions to the linear 


Cnm — 


dUn 


h P 

,m'^n 

(29) 

Cn{u). 

(30) 
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To obtain the stability matrix equation for the phase we Rrst derive the phase dynamical equation 
which corresponds to equation (5) for the modulus and which determines the phase variation 8(j>„ 
around the stationary solution 

exp exp +/( 5„_4 

^n'^n + l'^n + 2 exp ( T 

+ exp + 6(j>„+i)) 

+ (1 - e)A;„_ 2 Mn-iMn- 2 exp + (5(?i„_2)). (31) 


Here we have dropped the index 
of (31) gives 


zero at the stationary solution 
— = -T.{Dnm - C^^)8(j)m 


for convenience. Linearization 

(32) 


with 


Ct^ = 


dUr 


-Cn{u). 


(33) 


A similarity transformation of Dnm — C^^ with the matrix SnmUm leaves the eigenvalues unchanged 
and simplifies equation (32) to 


d 

dt 


n 


^(-^nm Cnm)^4^n 


(34) 


which only differs by a sign from equation (28). 

Both of the stability matrices —Dnm — C'nm and —Dnm + C'nm have their own set of eigen¬ 
values and eigenstates. As usual the eigenvalues, called cr, may be real or complex. The complex 
eigenvalues do not imply anything about the complexity of the perturbations around the static- 
velocity, which must be real. Rather a pair of complex conjugate eigenvalues represents a pair of 
perturbations which convert up and back into one another as time goes on and thereby produce 
an oscillatory response. 

To get numerical results, the eigenvalues of the non Hermitian Jacobians are calculated with 
the EISPACK package which is available via netlib. To achieve sufhcient numerical accuracy for 
the eigenvectors quadruple precision turns out to be necessary. 

When the real part of the eigenvalue is smaller than zero, the perturbation produced by the 
corresponding eigenstate decays. Conversely, when the real part is positive, the perturbation 
produces an exponentially growing response. For this reason, we focus our attention upon the 
eigenvalues with the largest real part. 

There are a large number, 2 • #, of eigenvalues. Figure 7 shows the pattern of eigenvalues 
for a typical situation (A = 2, # = 104, e = 0.3, v = 10“^ • 16“^® = 1.32349 • 10“^®). Notice 
how the eigenvalues arrange themselves into families. Each right eigenstate has a peak at some 
value of n (Fig. 8 b). Call this value m. Since the eigenvalues represent a relaxation rate, it is 
not surprising that for values of m toward the middle of the inertial range both the eigenvalues 
of the two branches in figure 7 have the magnitude a ~ as predicted by K41 scaling. The 

real family of eigenvalues continues into the dissipation range where a ~ since this is the 
magnitude of the dissipative relaxation rate. Only a few of the eigenvalues fail to fall into the 
pattern described here. These special eigenvalues have right eigenstates which tend to have a large 
weight at low shells. 

One phase disturbance has eigenvalue zero. This is reflective of a remaining phase ambiguity 
[13, 14] in the solution. Starting from any solution of our static equations one can obtain another 
one by the following process: All the u's with the same triality as M 5 are modified by multiplication 
by (with (f) real) while all the ones with the triality of ug are given the opposite phase. The 
new set of u„ is also a solution. This ambiguity in the solution is reflected in a zero eigenvalue, 
which will have no further importance for us. 


9 



For most values of A (the only exception being around A = 2.8) the largest real part occurs in 
the phase perturbation so that the phases show a greater tendency toward instability. As a result, 
the stability bounds shown in Rgure 1, which describe where the first eigenvalues go unstable, 
are phase instabilities. Each of these is also a pair of complex conjugate eigenvalues. Thus, the 
leading instabilities of the static solution are actually Hopf bifurcations in the phase variables. 
One such bifurcation gives periodic oscillatory behavior. According to the Ruelle Takens analysis, 
two such Hopf bifurcations can give chaotic behavior. Evidently, the Hopf bifurcations can lead to 
the GOY-model chaos. Indeed, this scenario has been demonstrated by Biferale et al. [12] for the 
GOY model. 

It is remarkable that the eigenvalue spectrums of the phase stability matrix C'nm — Dnm and 
of the (negative) magnitude stability matrix Cnm + Dnm are so different as the two matrices only 
differ by the dissipation contribution (29). This again is an indication for our main point of this 
paper: Dissipation effects can considerable change dynamical properties even in the inertial range. 

Luca Biferale [19] predicted that the GOY dynamics should turn unstable when the magnitude 
of the summands \U„\^(e— 1)“" in the second conserved (in the inviscous, unforced limit) quantity 
[^] Xln l^nP(£ ~ 1)“" increases with n. Plugging in K41 scaling, one obtains that the borderline 
of stability is 

\ = (35) 

We test this prediction by plotting the line of first and second phase instability in a log-log plot of 
A against 1 — e, cf. figure 9. Equation (35) catches the rough shape of the instability curves. This 
log-log plot also reveals that the stability curve crosses the curve of conserved helicity. It seems to 
oscillate around it and might asymptotically approach it. 

The largest real parts of the eigenvalues of the phase Jacobian are shown in figure 10, again 
for the standard parameters ly = 10“^, \ = 2, N = 22, but any larger N will lead to the same 
result, see above. In figure 11 we show how the eigenvalues a for the phase Jacobian matrix move 
within the complex cr-plane for increasing e. The above mentioned Hopf bifurcation is clearly seen. 
As mentioned in section 2, the solution ceases to exist around Cgnb ~ 0.4499. For the e-range 
]0, 0.4499] shown in figure 10, the modulus equations are linearly stable for these parameters. Yet 
this is not generally true: e.g. for lo = 10“^, A = 1.13, N = 104 (or larger), instability occurs 
around e^^^j Ri 0.11 in the phase equations and around e^^^j Ri 0.29 in the modulus equations. 

In figure 12 we present Cgtab = ^ttab ^ function of Igio for A = 2. Strong oscillations occur 
which mirror the oscillations of the with Ig discussed in section 3. Therefore, we expect a 
period of A'^, which indeed is the case. Our finding also suggests that the inertial range scaling 
exponents (p of the p-th velocity moments depend on the viscosity lo, as deviations 6(p = (p — pjS 
from the K41 value (p = p/3 set in smoothly at ef(aj(i^)- Such a dependence was already discussed 
in [4, 9]. 

The small wiggles on the curve in figure 12 are real. We only understand them in part. That 
every second oscillation should look the same can be understood along Biferale’s above mentioned 
argument [19] for the onset of instability: It will make a difference whether viscosity cuts in a shell 
with positive \U„\^(e — 1)“" (the summands of the second conserved quantity) or negative one, 
i.e., odd or even n. 

In figure 1 we have already presented and as a function of A for lo = 10“^. If turned by 
90 degrees, this figure resembles figure 4, the increasing fluctuations of and with increasing 

A. 

For larger A the fluctuations of become so large that stability and unstability change more 
than once with increasing e [14]. 

6 Summary and conclusions 

We set out in this paper to understand qualitatively the complex phase diagram of figure 1. We 
achieved that goal by explaining the observed oscillations in that figure in terms of a particular 
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and peculiar form of interaction between the integral scale and the VSR. 

Thus, in this paper, we have discussed a special coupling among the different turbulent scales. 
In this form of coupling, the integral scale produces an energy flux but does not fully determine 
its value. Instead the flux-value is set by the joint effect of the production and the dissipation 
mechanisms, i.e., though describing the ISR, it is set in the VSR. In the present calculations for 
the GOY model, the ISR-VSR coupling was made manifest by oscillations of the flux (and maybe 
also of inertial range scaling exponents (p) as a function of viscosity. These oscillations in turn 
result from dividing the wavevector space in Rnite shells and vanish only for diminishing shell 
distance. In shell models the oscillations are unavoidable, though we expect them to be much 
less pronounced (if noticeable at all) in the REWA calculations [11] as many more modes and 
couplings are present. In real turbulence there are no shells and no such oscillations. Nonetheless 
the more general lesson of this calculation might apply and there might well be a similar cooperative 
behavior between the VSR and the ISR, since both the corrections to scaling (if they exist) and 
the multifractal components (if they exist) are similarly undetermined by the inertial range as is 
the energy flux and the scaling corrections in the GOY model. All these quantities are set by the 
action of both the integral scale and the dissipative scale - working together. An independent hint 
for such a mechanism is experimentally found by Gastaing and collaborators [5] who clearly detect 
a Reynolds number dependence of ISR scaling exponents. 

Thus turbulent properties might well be a product of the detailed way the viscous range fits 
onto the ISR and turbulence might in end have a richer dependence upon viscous properties than 
contemplated in K41. If this finding can be confirmed, it will have prime importance for all models 
and simulations where the small scales are not fully resolved, as for example in hyperviscosity 
calculations or large eddy simulations. 
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Tables 



N = 3k + 2 

N = 3k + l 

II 

CO 

M4 



~ 1 

M 5 

~ 1 


-Vf 

Me 


~ 1 


^diss 

-P 

-P 


Ml 



- 

U2 

~ 1 


- 

M3 


~ 1 

- 

existence 

[ 0 ,oo[ 

[ 0 ,oo[ 

0 

0 

0 

0 

stability 

[ 0 ,oo[ 

[ 0 ,« 12 ] 

[ 0 ,f« 800] 


Table 1 

Scaling properties of the solutions for the three different cases N = ‘ik-\-2,N = ‘ik-\-l,N = ‘ik for 
large /. The hrst three lines follow from equations (24) and (25) and (23), the next one denotes 
the scaling of the total energy dissipation ediss- The scaling of ui, U 2 , and M 3 follows from the hrst 
three equations of (5). Note that in the last case N = 3k no solution exists. Correspondingly, the 
solution for the whole system of the N real equations for m„ ceases to exist for sufhciently large /. 
The approximate domain of existence and stability of the solution is given in the last two lines. 
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Figures 



Figure 1: Phase diagram which shows regions of existence and stability of the static solution. The 
right solid line is e^nb, left of which the static solution exists. Note that the minima and maxima of 
esni('^) roughly agree with those of J(A), cf. hgure 4. For completeness we mention that far right 
of this line there are (unmarked) regions where a stable solution of the GOY model again exists. 
The left solid line is (hrst phase instability), the middle thin solid line is (hrst magnitude 
instability). The calculation is done for v = 10“^ and N sufhciently large so that its value does not 
matter. The dotted line marks e = 1 — 1/A along which the helicity is conserved and along which 
dynamical scaling exponents C,q (beyond some threshold) are the same up to numerical accuracy 
[9]. The cross marks the standard case A = 2, e = 0.5. 
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Figure 2: Spectrum Un vs kn for e = 0.3, A = 2 for various different viscosities i/ = i/q = 10“^, 
V = vq ■ 16“^^, V = vq ■ 16“^^, V = vq ■ 16“^®, V = vq ■ 16“^®, left to right. Note that the detailed 
wiggle structure in the ISR, which now is identical for all cases, will be different when changing v 
by a factor different from (A'^)“™ where m is an integer. 
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Figure 3: (a) The energy flux J = for n in the ISR, showing strong oscillations with Igi^. The 
period is A'^. Here we chose N = 104, A = 2, e = 0.3. Note that the solution is unstable in some 
regions which we will discuss in detail in section 5. (b) The ratios (lower) and 

(upper) for k = 3, showing strong oscillations as J„. For all 3k in the ISR, these ratios are the 
same. 
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Figure 4: (a) The energy flux for n = 4, 11, 14, top to bottom. We chose N = 104, = 10“^, 

e = 0.3. The oscillations increase with increasing A. For large A viscous effects can already be felt 
for the large shells n = 11, 14. (b) Enlargement of the left part of the figure 4a in linear-log scale. 
The dashed line touching the maxima is 0.0016-f 0.0004IgA, according to (20). 
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Figure 5: Un vs / for n = 4, 6, 5 (top to bottom) for N = 22, c = 0.3, v = 10“^. Two regimes are 
identified: A small / regime with -^/-scaling and wiggles of period A®, and a large / regime with 
and scaling. 
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Figure 6: Energy spectrum Un vs kn in the large / regime: / m 10'^, e = 0.3, i/ = 10 = 22, 

_ 1/3 

A = 2. Also shown is the power law k„ ' . 
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Figure 7: Eigenvalues for the matrix equation (28) of the modulus for N = 104, e = 0.3, A = 2, 
v = 1.32349-10“^*^ in the complex plane. The complex pairs of eigenvalues (-Kcr, ±9 (t) are visualized 
as (Ig |iK(T|, lg(|9(T|), the real eigenvalues cr as (Ig \a\, 0). The largest eigenvalue cr Ri —10“® which is 
associated with the irregularity at M 3 is not seen in this plot. It is this eigenvalue that turns zero at 
fsni ~ 0.4499. Note that the irregularity in the pattern for (Ig \a\, 0) around 12 is not an numerical 
artefact but it is real and coincides with the end of the upper curve. Its physical meaning is the 
onset of the dissipation range VSR. 
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Figure 8: Eigenvector |jK(Mn)| vs. kn for magnitudes, N = 104, e = 0.3, \ = 2, v = 1.32349 • 10“^° 
for (a) the big eigenvalue a Ri —0.082 ± i0.035 (left and right eigenvectors) and (b) the smaller 
eigenvalue a Ri —2.0 • 10® ± il.O • 10^® (left and right eigenvectors). We did the calculation with 
quadrupole precision, the numerical noise level 10“^® can be recognized in both Rgures (a) and 
(b). 
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Figure 9: First and second phase (Hopf doublet) instabilities in a Ig-lg plot of A against 1 — e. 
Also shown is Biferale’s prediction [19], equation (35) (solid), and the line A = 1/(1 — e) along 
which “helicity” is the second conserved quantity [9] (dotted). The inset shows an enlargement 
of the small A range which shows that the instability line crosses the line of conserved “helicity”. 
The calculated points for the Rrst eigenvalue crossing zero are marked by dots; the line without 
dots denotes the second phase instability. 
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Figure 10: Largest eigenvalues of the phase Jacobian for the standard parameter case. The system 
turns unstable for e Ri 0.3699. A second mode becomes unstable at e Ri 0.3956. The capital letters 
refer to Rgure 11. 
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Figure 11: Movement of the eigenvalues a (of the phase Jacobian) in the complex plane for in¬ 
creasing e. We chose N = 22, = 10“^, A = 2. The Hopf bifurcation at e Ri 0.3699 is clearly 

identihed. The letters A, B, C denote the same eigenvalues as in Rgure 10. Eigenvalue D in that 
figure has much larger imaginary part 9(7 Ri 2 when crossing the imaginary axis and can thus also 
not be seen in this figure. Eigenvalue E in figure 10 (between zero and Ri 0.23) has zero imaginary 
part and can thus not be seen here. For this plot we chose Ae = 0.01. The lateral point density in 
each line of the plot is proportional to the inverse “speed” dMa/de in figure 10. The four arrows 
correspond to e = 0.1, 0.2, 0.3, and 0.4. 
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Figure 12: Dependence of Cgtab = ^ttab ^ standard parameter case A = 2. The shell 

number N is always chosen large enough so that we are in the inhnite shell number case. Note 
that these fluctuations are much smaller for smaller A, similar to those in the current J„, cf. Rgures 
3 and 4. They cannot account for the slight deviations between the Biferale prediction (35) and 
our numerical results in figure 9. 
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